
#include "Field.h"
#include "dbg.h"
#include "Cell.h"
#include "MemDebug.h"

template <int M>
void BuildCellData(
    const InputFile& file,
    std::vector<boost::shared_ptr<CellData<NData,M> > >& celldata, double& )
{
    const int n = file.getPos().size();
    celldata.resize(n);
    for(int i=0;i<n;++i) {
        celldata[i].reset(new CellData<NData,M>(file.getPos()[i]));
    }
}

template <int M>
void BuildCellData(
    const InputFile& file,
    std::vector<boost::shared_ptr<CellData<KData,M> > >& celldata, double& )
{
    Assert(file.getKappa().size() == file.getPos().size());
    Assert(file.getWeight().size() == 0 || 
           file.getWeight().size() == file.getPos().size());
    const int n = file.getPos().size();
    celldata.reserve(n);
    celldata.resize(0);
    for(int i=0;i<n;++i) {
        double w = file.getWeight().size() > 0 ? file.getWeight()[i] : 1.;
        if (w != 0.)
            celldata.push_back(boost::shared_ptr<CellData<KData,M> >(
                    new CellData<KData,M>(file.getPos()[i],file.getKappa()[i],w)));
    }
}

template <int M>
void BuildCellData(
    const InputFile& file,
    std::vector<boost::shared_ptr<CellData<EData,M> > >& celldata, double& vare)
{
    Assert(file.getShear().size() == file.getPos().size());
    Assert(file.getWeight().size() == 0 || 
           file.getWeight().size() == file.getPos().size());
    const int n = file.getPos().size();
    celldata.reserve(n);
    celldata.resize(0);
    double sumw = 0.;
    vare = 0.;
    for(int i=0;i<n;++i) {
        double w = file.getWeight().size() > 0 ? file.getWeight()[i] : 1.;
        if (w != 0.) {
            celldata.push_back(boost::shared_ptr<CellData<EData,M> >(
                    new CellData<EData,M>(file.getPos()[i],file.getShear()[i],w)));
            vare += w*w*norm(file.getShear()[i]);
            sumw += w;
        }
    }
    vare /= 2.*sumw*sumw/n; // 2 because we want var per component.
}


template <int DC>
Field<DC>::Field(const InputFile& file, const ConfigFile& params,
                 double minsep, double maxsep, double b) :
    _filename(file.getFileName()), _vare(0.)
{
    dbg<<"Building Field from file "<<file.getFileName()<<std::endl;

    // Read split_method as a string.
    std::string smstr = params.read<std::string>("split_method","mean");
    SplitMethod sm;
    if (smstr == "mean") sm = MEAN;
    else if (smstr == "median") sm = MEDIAN;
    else if (smstr == "middle") sm = MIDDLE;
    else {
        myerror("Invalid split_method",smstr.c_str());
        sm = MEAN; // To prevent compiler warning
    }
    xdbg<<"Splitting using "<<smstr<<std::endl;

    // We don't bother accumulating the mean information for Cells that would be 
    // too large or too small.

    // The minimum size cell that will be useful is one where two cells that just barely
    // don't split have (d + s1 + s2) = minsep
    // i.e. d = minsep - 2s  and 2s = bd
    //      d = minsep - bd
    //      d = minsep / (1+b)
    //      s = b/2 * minsep / (1+b)
    double minsize = 0.5 * minsep * b / (1.+b);
    double minsizesq = minsize * minsize;
    xdbg<<"minsizesq = "<<minsizesq<<std::endl;

    // The maximum size cell that will be useful is one where a cell of size s will
    // be split at the maximum separation even if the other size = 0.
    // i.e. s = b * maxsep
    double maxsize = maxsep * b;
    double maxsizesq = maxsize * maxsize;
    xdbg<<"maxsizesq = "<<maxsizesq<<std::endl;

    // Build the right kind of Cell according to whether the InputFile uses RA,Dec or x,y.
    if (file.useRaDec()) {
        _metric = Sphere;
        std::vector<boost::shared_ptr<CellData<DC,Sphere> > > celldata;
        BuildCellData(file,celldata,_vare);
        xdbg<<"Built celldata with "<<celldata.size()<<" entries\n";
        if (DC == EData)
            dbg<<"vare = "<<_vare<<": sig_sn (per component) = "<<sqrt(_vare)<<std::endl;
        buildCells(celldata,_cells_sphere,minsizesq,maxsizesq,sm,0,celldata.size());
    } else {
        _metric = Flat;
        std::vector<boost::shared_ptr<CellData<DC,Flat> > > celldata;
        BuildCellData(file,celldata,_vare);
        xdbg<<"Built celldata with "<<celldata.size()<<" entries\n";
        if (DC == EData)
            dbg<<"vare = "<<_vare<<": sig_sn (per component) = "<<sqrt(_vare)<<std::endl;
        buildCells(celldata,_cells_flat,minsizesq,maxsizesq,sm,0,celldata.size());
    }
}

template <int DC> template <int M>
void Field<DC>::buildCells(
    std::vector<boost::shared_ptr<CellData<DC,M> > >& vdata,
    std::vector<boost::shared_ptr<Cell<DC,M> > >& cells,
    double minsizesq, double maxsizesq,
    SplitMethod sm, size_t start, size_t end)
{
    //xxdbg<<"Start buildCells: start,end = "<<start<<','<<end<<std::endl;
    //xxdbg<<"minsizesq = "<<minsizesq<<", maxsizesq = "<<maxsizesq<<std::endl;
    // The structure of this is very similar to the Cell constructor.
    // The difference is that here we only construct a new Cell (and do the corresponding
    // calculation of the averages) if the size is small enough.  At that point, the 
    // rest of the construction is passed onto the Cell class.
    boost::shared_ptr<CellData<DC,M> > ave;
    double sizesq;
    if (end-start == 1) {
        //xxdbg<<"Only 1 CellData entry: size = 0\n";
        ave = vdata[start];
        sizesq = 0.;
    } else {
        ave.reset(new CellData<DC,M>(vdata,start,end));
        sizesq = CalculateSizeSq(ave->pos,vdata,start,end);
        //xxdbg<<"size = "<<sqrt(sizesq)<<std::endl;
    }

    if (sizesq <= maxsizesq) {
        //xxdbg<<"Small enough.  Make a cell.\n";
        if (end-start > 1) ave->finishAverages(vdata,start,end);
        boost::shared_ptr<Cell<DC,M> > p(new Cell<DC,M>(ave,sizesq,vdata,minsizesq,sm,start,end));
        cells.push_back(p);
    } else {
        size_t mid = SplitData(vdata,sm,start,end,ave->pos);
        //xxdbg<<"Too big.  Recurse with mid = "<<mid<<std::endl;
        buildCells(vdata,cells,minsizesq,maxsizesq,sm,start,mid);
        buildCells(vdata,cells,minsizesq,maxsizesq,sm,mid,end);
    }
}

template class Field<NData>;
//template class Field<KData>;
template class Field<EData>;
